---
title: "Replication Material for Politics of Welfare Exclusion: Open or Concealed Welfare Chauvinism and Public Support for Social Assistance Reform"
author: "Mikkel Haderup Larsen"
date: "`r Sys.Date()`"
output:
  html_document:
    df_print: paged
    toc: yes
    toc_depth: '2'
  pdf_document:
    toc: yes
    toc_depth: 2
  word_document:
    toc: yes
    toc_depth: '2'
geometry: margin=2cm
fontsize: 11pt
---

\newpage
# Introduction - Read Me
This online supplement contains results of all analyses and the R code that produces these results based on the raw data that is also part of the online supplement. If you are interested in the results and maybe also in skimming the R code that generates them, you can simply go through this document.

However, if you wish to replicate the results themselves and maybe alter the R code, you will need to work with the "Chauv_policy.rmd" file. To use the Rmd replication file, you need a working version of R (see detailed session info under "Setup" below), RStudio, RMarkdown and Latex. To install RMarkdown and Latex, use the following code after you have installed R and RStudio:

```{r eval = FALSE}
install.packages(c('rmarkdown', 'tinytex'), dep = TRUE)
tinytex::install_tinytex()  # install TinyTeX
```

All mentioned software packages are open source, free to use and available for Windows, Mac Os, and Unix operating systems. Once you have installed these, please open RStudio and establish an RStudio project in the unzipped file directory that contains the ".Rmd" and ".dta" files. For further information on Rstudio projects and how to set them up, please refer to [R for Data Science, Chapter 8](https://r4ds.had.co.nz/workflow-projects.html). 

After having set up the RStudio project, you will find all necessary files under the "Files" tab in Rstudio. Open the file "Chauv_policy.Rmd" by clicking on it under the "Files" tab. The file should open in the top left panel. Rstudio will automatically notify you about all user-written packages that the analysis relies on and that you have not yet downloaded and installed. Please click on "install". To replicate the whole analysis, simply click on "Knit", just below the file tab of the top left panel showing the Rmd file. If you want to alter the R code, you can work on the single R-snippets. For more information on how to edit Rmd files and R snippets see [R for Data Science, Chapter 27](https://r4ds.had.co.nz/r-markdown.html).

\newpage
```{r knitr, include = FALSE}
# Knitr options for this document
options(
  tinytex.verbose = TRUE, 
  digits = 3)

knitr::opts_chunk$set(
  warning = FALSE, message = FALSE, cache = FALSE,
  collapse = TRUE, fig.retina = 3, 
  fig.align='center', error = TRUE, digits = 2)
```

In order to perfectly reproduce the results, you might need to ensure that you use the same versions of R and of the R-packages I used. The session info at the bottom of this code chunk provides all necessary information.

\newpage
# Load libraries
```{r Library}
# Libraries
pacman::p_load(haven,
               tidyverse,
               lubridate,
               stringr,
               survey,
               effects,
               estimatr,
               cregg,
               cjbart,
               marginaleffects,
               cowplot,
               ggplot2,
               ggpubr,
               texreg,
               gtsummary,
               kableExtra) 

# Show session info on used packages and versions
sessionInfo()
```

\newpage
# Load, clean, and tidy data
```{r Data}
Benefits_data <- read_dta("INSERT PATH HERE") 

Benefits_data <- Benefits_data %>%
mutate(
ID = CBC_Id %>% as.numeric(),
Choice = case_when( ## Binary outcome: was the social assistance reform proposal chosen?
CBC_Random1 == 1 & concept == 1 & task == 1 |
CBC_Random2 == 1 & concept == 1 & task == 2 |
CBC_Random3 == 1 & concept == 1 & task == 3 |
CBC_Random4 == 1 & concept == 1 & task == 4 |
CBC_Random1 == 2 & concept == 2 & task == 1 |
CBC_Random2 == 2 & concept == 2 & task == 2 |
CBC_Random3 == 2 & concept == 2 & task == 3 |
CBC_Random4 == 2 & concept == 2 & task == 4 ~ 1,
CBC_Random1 == 1 & concept == 2 & task == 1 |
CBC_Random2 == 1 & concept == 2 & task == 2 |
CBC_Random3 == 1 & concept == 2 & task == 3 |
CBC_Random4 == 1 & concept == 2 & task == 4 |
CBC_Random1 == 2 & concept == 1 & task == 1 |
CBC_Random2 == 2 & concept == 1 & task == 2 |
CBC_Random3 == 2 & concept == 1 & task == 3 |
CBC_Random4 == 2 & concept == 1 & task == 4 ~ 0),
Proposer = case_when( ## The treatments
    attr01forslagstiller == 1 ~ "Party from red bloc",
    attr01forslagstiller == 2 ~ "Party from blue bloc") %>% factor() %>% fct_relevel("Party from blue bloc",
                                         "Party from red bloc"),
Work = case_when(
    Attr02Kravomudførtarbejdeindenfo == 1 ~ "No work requirements",
    Attr02Kravomudførtarbejdeindenfo == 2 ~ "225 hours",
    Attr02Kravomudførtarbejdeindenfo == 3 ~ "225 hours both partners in a union",
    Attr02Kravomudførtarbejdeindenfo == 4 ~ "300 hours both partners in a union") %>% factor() %>%
  fct_relevel("No work requirements",
              "225 hours",
              "225 hours both partners in a union",
              "300 hours both partners in a union"),
Residence = case_when(
    Attr03Bopælskrav == 1 ~ "Residency in DK",
    Attr03Bopælskrav == 2 ~ "Lived in DK 3 out of 4 years",
    Attr03Bopælskrav == 3 ~ "Lived in DK 9 out of 10 years",
    Attr03Bopælskrav == 4 ~ "Lived in DK 10 out of 11 years")  %>% factor() %>%
    fct_relevel("Residency in DK",
                "Lived in DK 3 out of 4 years",
                "Lived in DK 9 out of 10 years",
                "Lived in DK 10 out of 11 years"),
Citizenship = case_when(
    Attr04Kontanthjælpsatsforpersone == 1 ~ "Non-citizens receive 50% less",
    Attr04Kontanthjælpsatsforpersone == 2 ~ "Non-citizens receive 25% less",
    Attr04Kontanthjælpsatsforpersone == 3 ~ "Non-citizens receive the same") %>% factor() %>%
  fct_relevel("Non-citizens receive the same",
              "Non-citizens receive 25% less",
              "Non-citizens receive 50% less"),
Handicap = case_when(
    Attr05Kontanthjælpssatsforperson == 1 ~ "Disabled receive 50% more",
    Attr05Kontanthjælpssatsforperson == 2 ~ "Disabled receive 25% more",
    Attr05Kontanthjælpssatsforperson == 3 ~ "Disabled receive the same") %>% factor() %>%
  fct_relevel("Disabled receive the same",
              "Disabled receive 25% more",
              "Disabled receive 50% more"),
Parents = case_when(
    Attr06Kontanthjælpsatsforforældr == 1 ~ "Parents receive 50% more",
    Attr06Kontanthjælpsatsforforældr == 2 ~ "Parents receive 25% more",
    Attr06Kontanthjælpsatsforforældr == 3 ~ "Parents receive the same") %>% factor() %>%
  fct_relevel("Parents receive the same",
              "Parents receive 25% more",
              "Parents receive 50% more"),
Sanctions = case_when(
    Attr07Sanktionvedudeblivelsefrae == 1 ~ "No sanctions",
    Attr07Sanktionvedudeblivelsefrae == 2 ~ "Receive 10% less the following month",
    Attr07Sanktionvedudeblivelsefrae == 3 ~ "Receive 20% less the following month") %>% factor() %>%
  fct_relevel("No sanctions",
              "Receive 10% less the following month",
              "Receive 20% less the following month"),
Joboffer = case_when(
    attr08kravomacceptafjobtilbudder == 1 ~ "Required to accept joboffer that matches skills",
    attr08kravomacceptafjobtilbudder == 2 ~ "Not required to accept joboffer that matches skills") %>% factor() %>%
  fct_relevel("Not required to accept joboffer that matches skills",
              "Required to accept joboffer that matches skills"),
Jobtraining = case_when(
    Attr09Kravomdeltagelseijobtrænin == 1 ~ "Required to participate in job training",
    Attr09Kravomdeltagelseijobtrænin == 2 ~ "Not required to participate in job training") %>% factor() %>%
  fct_relevel("Not required to participate in job training",
              "Required to participate in job training"),
Costs = case_when(
    Attr10Årligskattestigningsomfølg == 1 ~ "Lower personal costs than previous reform",
    Attr10Årligskattestigningsomfølg == 2 ~ "Same personal costs as previous reform",
    Attr10Årligskattestigningsomfølg == 3 ~ "Higher personal costs than previous reform") %>% factor() %>%
  fct_relevel("Same personal costs as previous reform",
              "Lower personal costs than previous reform",
              "Higher personal costs than previous reform"),
  Gender = case_when( ### The controls
      bagg2 == 1 ~ "Male",
      bagg2 == 2 ~ "Female") %>% factor() %>% fct_relevel("Male",
                                                          "Female"),
  Age = case_when(
    bagg1_a == 2 ~ 24,
    bagg1_a == 9 ~ 70,
    bagg1_o1 < 18 ~ NA_real_,
    TRUE ~ zap_labels(bagg1_o1)),
  University = case_when(
    bagg5 < 7 ~ "No university degree",
    bagg5 > 6 & bagg5 < 9  ~ "University degree") %>% factor() %>%
    fct_relevel("No university degree",
                "University degree"),
  Copenhagen = case_when(
      Region != 4 ~ "Other area",
      Region == 4 ~ "Greater Copenhagen area") %>% factor() %>%
    fct_relevel("Other area",
                "Greater Copenhagen area"),
  White_collar = case_when(
    bagg6 != 3 ~ "Not white collar",
    bagg6 == 3 ~ "White collar") %>% factor() %>%
    fct_relevel("Not white collar",
                "White collar"),
  Blue_collar = case_when(
    bagg6 > 2 ~ "Not blue collar",
    bagg6 < 3 ~ "Blue collar") %>% factor() %>%
    fct_relevel("Not blue collar",
                "Blue collar"),
  Retired = case_when(
    bagg6 == 8 ~ "Retired",
    bagg6 != 8 ~ "Not retired")  %>% factor() %>%
    fct_relevel("Not retired",
                "Retired"),
  Income = case_when(
    bagg7 < 8  ~ "Reported annual income less than 400,000 DKK",
    bagg7 > 7 & bagg7 < 16 ~ "Reported annual income more than 400,000 DKK") %>% factor() %>%
    fct_relevel("Reported annual income less than 400,000 DKK",
                "Reported annual income more than 400,000 DKK"),
  Device = case_when(
    sys_OperatingSystem == "Android 10" |
    sys_OperatingSystem == "Android 11" |
    sys_OperatingSystem == "Android 12" |
    sys_OperatingSystem == "Android 5.0" |
    sys_OperatingSystem == "Android 5.1.1" |
    sys_OperatingSystem == "Android 7.0" |
    sys_OperatingSystem == "Android 7.1.1" |
    sys_OperatingSystem == "Android 8.0.0" |
    sys_OperatingSystem == "Android 8.1.0" |
    sys_OperatingSystem == "Android 9" |
    sys_OperatingSystem == "iPad" |
    sys_OperatingSystem == "iPhone" ~ "Phone/Tablet",
    sys_OperatingSystem == "Chrome OS" |
    sys_OperatingSystem == "Linux" |
    sys_OperatingSystem == "Mac OS X 10.10" |
    sys_OperatingSystem == "Mac OS X 10.10.5" |
    sys_OperatingSystem == "Mac OS X 10.12.6" |
    sys_OperatingSystem == "Mac OS X 10.13" |
    sys_OperatingSystem == "Mac OS X 10.13.6" |
    sys_OperatingSystem == "Mac OS X 10.14.6" |
    sys_OperatingSystem == "Mac OS X 10.15" |
    sys_OperatingSystem == "Mac OS X 10.15.2" |
    sys_OperatingSystem == "Mac OS X 10.15.5" |
    sys_OperatingSystem == "Mac OS X 10.15.6" |
    sys_OperatingSystem == "Mac OS X 10.15.7" |
    sys_OperatingSystem == "Windows 7" |
    sys_OperatingSystem == "Windows 8" |
    sys_OperatingSystem == "Windows 8.1" |
    sys_OperatingSystem == "Windows NT 10.0" ~ "PC/Mac") %>% factor() %>% fct_relevel("Phone/Tablet", "PC/Mac"),
  Fear_unemployment = case_when(
    Q1 == 1 | Q1 == 2 | bagg6 == 5 | bagg6 == 6 | bagg6 == 8 | bagg6 == 9 ~ "No",
    Q1 == 3 | Q1 == 4 | bagg6 == 7 ~ "Yes") %>% factor() %>% fct_relevel("No","Yes"),
  Bloc = case_when(
    Q2 == 1 | Q2 == 2 | Q2 == 5 | Q2 == 6 | Q2 == 11 | Q2 == 13 | Q2 == 14 ~ "Voted red",
    Q2 == 3 | Q2 == 4 | Q2 == 7 | Q2 == 8 | Q2 == 9 | Q2 == 10 | Q2 == 12 ~ "Voted blue") %>% factor() %>% fct_relevel ("Voted red", "Voted blue"),
  News_welfare = case_when(
    Q3 == 1 | Q3 == 2 |Q3 == 3 ~ "Reads news about welfare state at least once a week",
    Q3 > 3 ~ "Reads news about welfare state less than once a week") %>% factor() %>% fct_relevel("Reads news about welfare state less than once a week",
                                                            "Reads news about welfare state at least once a week"),
  SurveyDate = str_extract(ActualSurveyEndTime, "\\d{4}-\\d{2}-\\d{2}"),
  Day = case_when(
      SurveyDate == "2022-05-09" ~ 1,
      SurveyDate == "2022-05-10" ~ 2,
      SurveyDate == "2022-05-11" ~ 3,
      SurveyDate == "2022-05-12" ~ 4,
      SurveyDate == "2022-05-13" ~ 5,
      SurveyDate == "2022-05-14" ~ 6,
      SurveyDate == "2022-05-15" ~ 7,
      SurveyDate == "2022-05-16" ~ 8,
      SurveyDate == "2022-05-17" ~ 9,
      SurveyDate == "2022-05-18" ~ 10,
      SurveyDate == "2022-05-19" ~ 11) %>% as.numeric(),
  Start_time = ActualSurveyStartTime,
  End_time = ActualSurveyEndTime,
  Responsetime = as.numeric(difftime(End_time, Start_time, units = "mins")),
  Profile = case_when(
      concept == 1 ~ "Left",
      concept == 2 ~ "Right") %>% factor() %>%
    fct_relevel("Left",
                "Right"),
    Task = case_when(
      task == 1 ~ "First choice task",
      task == 2 ~ "Second choice task",
      task == 3 ~ "Third choice task",
      task == 4 ~ "Fourth choice task") %>% factor() %>%
    fct_relevel("First choice task",
                "Second choice task",
                "Third choice task",
                "Fourth choice task"),
    Weight = wgt,
Row_order = str_replace_all(qOrder, "[[:punct:]]", ",") %>% as.character())%>% select(ID, Task, Profile, task, Weight, Choice, Work, Residence, Citizenship, Proposer, Costs, Handicap, Parents, Sanctions, Joboffer, Jobtraining, Device, University, Age, Gender, Income, Blue_collar, White_collar, Retired, Copenhagen, Day, SurveyDate, Start_time, End_time, Responsetime, News_welfare, Fear_unemployment, Bloc, Row_order)

Benefits_data <- (subset(Benefits_data, ID != 530))

Benefits_analysis <- Benefits_data %>% select(ID, Task, Profile, task, Choice, Weight, Work, Residence, Citizenship, Proposer, Costs, Handicap, Parents, Sanctions, Joboffer, Jobtraining, Weight, Row_order) %>% drop_na()

Benefits_analysis <- Benefits_analysis %>%
  separate(Row_order, into = c("Row_x", "Row_1","Row_2","Row_3","Row_4","Row_5","Row_6","Row_7","Row_8","Row_9", "Row_10", "Row_y"), sep=",", extra = "merge") %>% select(-Row_x, -Row_y)

Benefits_analysis <- Benefits_analysis %>%
mutate(
  Residence_row = case_when(
    Row_1 == 3 ~ "1",
    Row_2 == 3 ~ "2",
    Row_3 == 3 ~ "3",
    Row_4 == 3 ~ "4",
    Row_5 == 3 ~ "5",
    Row_6 == 3 ~ "6",
    Row_7 == 3 ~ "7",
    Row_8 == 3 ~ "8",
    Row_9 == 3 ~ "9",
    Row_10 == 3 ~ "10") %>% factor() %>% fct_relevel("1","2", "3", "4", "5", "6", "7", "8", "9", "10"),
   Citizenship_row = case_when(
    Row_1 == 4 ~ "1",
    Row_2 == 4 ~ "2",
    Row_3 == 4 ~ "3",
    Row_4 == 4 ~ "4",
    Row_5 == 4 ~ "5",
    Row_6 == 4 ~ "6",
    Row_7 == 4 ~ "7",
    Row_8 == 4 ~ "8",
    Row_9 == 4 ~ "9",
    Row_10 == 4 ~ "10") %>% factor() %>% fct_relevel("1","2", "3", "4", "5", "6", "7", "8", "9", "10"))

Benefits_subgroup <- Benefits_data %>%
  mutate(
        Observation = row_number()) %>%
  select(Observation, Choice, Residence, Citizenship,Age, Gender, Copenhagen, University, Income, Blue_collar, White_collar, Retired, News_welfare, Fear_unemployment, Bloc, Day, Responsetime, Device) %>% drop_na()

Benefits_reduced <- Benefits_data %>%
  mutate(
        Observation = row_number()) %>%
  select(Observation,ID, Weight, Choice, Residence, Citizenship, Proposer, Costs, Handicap, Parents, Sanctions, Joboffer, Jobtraining, Work, Age, Gender, Copenhagen, University, Income, Blue_collar, White_collar, Retired, News_welfare, Fear_unemployment, Bloc, Day, Responsetime, Device) %>% drop_na()
```

\newpage
# Main analysis
## Regression tables (Appendix D)
```{r Tables, results = FALSE}
#--------------#
# H1-H2. AMCEs #
#--------------#
Main_mod <- lm_robust(
  data = Benefits_analysis, 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod)

Main_mod_reduced <- lm_robust(
  data = Benefits_reduced, 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod_reduced)

htmlreg(
  list(Main_mod, Main_mod_reduced),
  custom.header = list("Main analysis (Full sample)" = 1, "Main analysis (Reduced sample)" =2),
  include.ci = TRUE, include.rmse = FALSE, digits = 3,
  doctype = FALSE, stars = c(0.001, 0.01, 0.05, 0.1), symbol = "+", 
  custom.coef.names = c("(Intercept)", "Must have lived in DK 3 out of 4 years", "Must have lived in DK 9 out of 10 years", "Must have lived in DK 10 out of 11 years", "Non-citizens receive 25% less", "Non-citizens receive 50% less", "Proposed by party from red bloc", "Lower costs than status quo", "Higher costs than status quo", "Disabled receive 25% more", "Disabled receive 50% more", "Parents receive 25% more", "Parents receive 50% more", "No-show at employment agency meeting: 10% sanction for 1 month", "No-show at employment agency meeting: 20% sanction for 1 month", "Must accept job offer that matches skills", "Must attend job training", "Must have worked minimum 225 hours last 12 months","Must have worked minimum 225 hours last 12 months. Rule also applies to spouse", "Must have worked minimum 300 hours last 12 months. Rule also applies to spouse"),
file = "Table_D1.doc", star.symbol = "\\*")
```

\newpage
## Figure 1: Average marginal component effects

\newpage
Figure 1 displays the AMCEs of the manipulated policy attributes. The figure is identical to Figure 1 in the main article.
```{r Figure_1A, results='asis'}
#---------------------------
# H1-H2. "Average effects": AMCE
#---------------------------
## Generate the plotdata
plotdata1 <- tidy(Main_mod) %>% 
  filter(term != "(Intercept)") %>%
  mutate(
    min90 = (estimate - qt(0.95, df) * std.error) * 100,
    min95 = (estimate - qt(0.975, df) * std.error) * 100,
    max90 = (estimate + qt(0.95, df) * std.error) * 100,
    max95 = (estimate + qt(0.975, df) * std.error) * 100,
    Attribute = case_when(
      grepl("Residence", term) ~ "Residence",
      grepl("Non-citizens", term) ~ "Non-citizens",
      grepl("Party", term) ~ "Proposer",
      grepl("costs", term) ~ "Tax costs",
      grepl("Disabled", term) ~ "Disability",
      grepl("Parents", term) ~ "Parents",
      grepl("following", term) ~ "Sanctions",
      grepl("accept", term) ~ "Job offer",
      grepl("training", term) ~ "Job training",
      grepl("Work", term) ~ "Prior work") %>% factor() %>% fct_relevel("Residence",
                                               "Non-citizens",
                                               "Disability",
                                               "Parents",
                                               "Proposer",
                                               "Tax costs",
                                               "Sanctions",
                                               "Job offer",
                                               "Job training",
                                               "Prior work"),
    label_up = round(estimate, digits = 3) * 100,
    estimate = estimate * 100) %>% drop_na()

plotdata1 <- plotdata1 %>% add_row(term=c("Lives in DK", 
                                          "Non-citizens receive the same",
                                          "Disabled receive the same",
                                          "Parents receive the same",
                                          "Proposer is from blue bloc",
                                          "No change in annual taxes",
                                          "No sanctions",
                                          "No job offer requirements",
                                          "No job training requirements",
                                          "No work requirements"), 
                                   estimate=c(0,0,0,0,0,0,0,0,0,0), 
                                   Attribute=c("Residence",
                                               "Non-citizens",
                                               "Disability",
                                               "Parents",
                                               "Proposer",
                                               "Tax costs",
                                               "Sanctions",
                                               "Job offer",
                                               "Job training",
                                               "Prior work"))

plotdata1 <- plotdata1 %>%
  mutate(
    term = case_when(
      term == "No work requirements" ~ "No work requirements",
      term == "Work225 hours" ~ "225 h last 12 months",
      term == "Work225 hours both partners in a union" ~ "225 h last 12 months incl. spouse",
      term == "Work300 hours both partners in a union" ~ "300 h last 12 months incl. spouse",
      term == "Lives in DK" ~ "Lives in DK",
      term == "ResidenceLived in DK 3 out of 4 years" ~ "Lived in DK 3 of 4 years",
      term == "ResidenceLived in DK 9 out of 10 years" ~ "Lived in DK 9 of 10 years",
      term == "ResidenceLived in DK 10 out of 11 years" ~ "Lived in DK 10 of 11 years",
      term == "Non-citizens receive the same" ~ "Non-citizens receive the same",
      term == "CitizenshipNon-citizens receive 25% less" ~ "Non-citizens receive 25% less",
      term == "CitizenshipNon-citizens receive 50% less" ~ "Non-citizens receive 50% less",
      term == "Proposer is from blue bloc" ~ "Right-wing party",
      term == "ProposerParty from red bloc" ~ "Left-wing party",
      term == "CostsLower personal costs than previous reform" ~ "Pay 0.5% less in annual taxes",
      term == "No change in annual taxes" ~ "No change in annual taxes",
      term == "CostsHigher personal costs than previous reform" ~ "Pay 0.5% more in annual taxes",
      term == "Disabled receive the same" ~ "Disabled receive the same",
      term == "HandicapDisabled receive 25% more" ~ "Disabled receive 25% more",
      term == "HandicapDisabled receive 50% more" ~ "Disabled receive 50% more",
      term == "Parents receive the same" ~ "Parents receive the same",
      term == "ParentsParents receive 25% more" ~ "Parents receive 25% more",
      term == "ParentsParents receive 50% more" ~ "Parents receive 50% more",
      term == "No sanctions" ~ "No sanctions",
      term == "SanctionsReceive 10% less the following month" ~ "Receive 10% less for 1 month",
      term == "SanctionsReceive 20% less the following month" ~ "Receive 20% less for 1 month",
      term == "No job offer requirements" ~ "No job offer requirements",
      term == "JobofferRequired to accept joboffer that matches skills" ~ "Must accept job offer",
      term ==  "No job training requirements" ~  "No job training requirements",
      term == "JobtrainingRequired to participate in job training" ~ "Must attend job training") %>% factor() 
                                                                              %>% fct_relevel("300 h last 12 months incl. spouse",
                                                                                              "225 h last 12 months incl. spouse",
                                                                                              "225 h last 12 months",
                                                                                              "No work requirements",
                                                                                              "Must attend job training",
                                                                                              "No job training requirements",
                                                                                              "Must accept job offer",
                                                                                              "No job offer requirements",
                                                                                              "Receive 20% less for 1 month",
                                                                                              "Receive 10% less for 1 month",
                                                                                              "No sanctions",
                                                                                              "Left-wing party",
                                                                                              "Right-wing party",
                                                                                              "Pay 0.5% more in annual taxes",
                                                                                              "No change in annual taxes",
                                                                                              "Pay 0.5% less in annual taxes",
                                                                                              "Parents receive 50% more",
                                                                                              "Parents receive 25% more",
                                                                                              "Parents receive the same",                                                                                              
                                                                                              "Disabled receive 50% more",
                                                                                              "Disabled receive 25% more",
                                                                                              "Disabled receive the same",                                                                                               
                                                                                              "Non-citizens receive 50% less",
                                                                                              "Non-citizens receive 25% less",
                                                                                              "Non-citizens receive the same",
                                                                                              "Lived in DK 10 of 11 years",
                                                                                              "Lived in DK 9 of 10 years",
                                                                                              "Lived in DK 3 of 4 years",
                                                                                              "Lives in DK"))

plotdata1 <- plotdata1 %>%
  mutate(
    Attribute = Attribute %>% factor() %>% fct_relevel("Residence",
                                               "Non-citizens",
                                               "Disability",
                                               "Parents",
                                               "Proposer",
                                               "Tax costs",
                                               "Sanctions",
                                               "Job offer",
                                               "Job training",
                                               "Prior work"))


## Plot   
Figure_1 <- ggplot(data = plotdata1, 
       aes(y = estimate, x = term, 
           ymin = min90, ymax = max90)) +
  geom_hline(yintercept = 0, color = "#901A1E", lty = "dashed") +
  geom_linerange(aes(ymin = min95, ymax = max95), color = "#808080") +
  geom_pointrange(size = 0.5) +
  geom_text(aes(label = label_up), vjust = -1.2, hjust = -0.5, size = 3) +
scale_y_continuous(breaks = c(-10,-7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10),
labels = c("-10", "-7.5", "-5", "-2.5", "0", "2.5", "5", "7.5", "10")) +
coord_flip() +
  labs(x = "", 
       y = "Average Marginal Component Effect") + 
  facet_grid(Attribute~., scales = "free_y", switch = "y") +
  theme_classic2() 

## Save Figure 1
ggsave("Figure_1.jpg", plot = last_plot(), width = 8, height = 10, dpi = 600)
```

\newpage
```{r Figure_1B, results='asis'}
Main_mod <- lm_robust(
  data = Benefits_analysis, 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod)

Main_mod_importance <- cj(Benefits_analysis, Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, estimate = "amce",
                  id = ~ ID) %>% 
  group_by(feature) %>%
  summarise(max_estimate = max(estimate),
            min_estimate = min(estimate),
            range = max_estimate-min_estimate) %>%
  ungroup %>%
  mutate(
    sum_range = sum(range),
    importance = (range / sum_range) * 100,
    sum_importance = sum(importance)) %>% 
  mutate(
    wc = case_when(
      feature %in% c("Citizenship", "Residence") ~ "Welfare Chauvinism",
      feature %in% c("Proposer", "Costs", "Handicap", "Parents", "Sanctions", "Joboffer", "Jobtraining", "Work") ~ "Other attributes") %>% factor(),
    Attribute = case_when(
      feature == "Residence" ~ "Length-of-residence",
      feature == "Citizenship" ~ "Citizenship status",
      feature == "Proposer" ~ "Proposer",
      feature == "Costs" ~ "Tax costs",
      feature == "Handicap" ~ "Disability",
      feature == "Parents" ~ "At-home kids",
      feature == "Sanctions" ~ "Sanctions",
      feature == "Joboffer" ~ "Job offer",
      feature == "Jobtraining" ~ "Job training",
      feature == "Work" ~ "Prior work") %>% factor() %>% fct_relevel("Length-of-residence",
                                               "Citizenship status",
                                               "Disability",
                                               "At-home kids",
                                               "Proposer",
                                               "Tax costs",
                                               "Sanctions",
                                               "Job offer",
                                               "Job training",
                                               "Prior work"))
summary(Main_mod_importance)

Figure_importance <- ggplot(Main_mod_importance, aes(fct_reorder(Attribute, importance, .desc = F), importance, fill = wc)) + 
  geom_col() +
  coord_flip() +
  labs(x = NULL, y = "Importance in %") +
  scale_fill_manual(
    values = c(
      "Welfare Chauvinism" = "#000000",
      "Other attributes" = "#808080")) +
  scale_y_continuous(limits = c(0,20), breaks = c(0,5,10,15,20)) +
  theme_classic2() +
  theme(legend.position = "bottom",
        legend.title = element_blank())

ggsave("Figure_importance.jpg", plot = last_plot(), width = 6, height = 4, dpi = 400)

```

\newpage
## Figure 2: Sub-group effects 
Figure 2 displays the results of the sub-group analysis using Bayesian Additive Regression Trees. The figure is identical to Figure 2 in the main article.
```{r Figure_2, results='asis'}

BART_model <- cjbart(data = Benefits_subgroup,
                   Y = "Choice", 
                   id = "Observation",
                   cores = 4, 
                   seed = 89)

IMCE_model <- IMCE(Benefits_subgroup,
                   BART_model,
                   attribs = c("Residence", "Citizenship"),
                   ref_levels = c("Residency in DK",
                                  "Non-citizens receive the same"),
                       cores = 4) 

# VIMP analysis
#IMCE_model[sapply(IMCE_model, is.character)] <- lapply(IMCE_model[sapply(IMCE_model, is.character)],as.factor)
Hetereogeneity_results <- het_vimp(IMCE_model,
                                   levels = c("Non-citizens receive 25% less",
                                              "Non-citizens receive 50% less",
                                              "Lived in DK 3 out of 4 years",
                                              "Lived in DK 9 out of 10 years",
                                              "Lived in DK 10 out of 11 years"),
                                   covars = c("Age",
                                              "Gender",
                                              "Copenhagen",
                                              "University",
                                              "Income",
                                              "Blue_collar",
                                              "White_collar",
                                              "Retired",
                                              "News_welfare",
                                              "Fear_unemployment",
                                              "Bloc",
                                              "Day", 
                                              "Responsetime",
                                              "Device")) # Very computationally demanding

# Tidy labels
Covariate_lookup <- data.frame(Covariate_name = c("Age",
                                              "Gender",
                                              "Copenhagen",
                                              "University",
                                              "Income",
                                              "Blue_collar",
                                              "White_collar",
                                              "Retired",
                                              "News_welfare",
                                              "Fear_unemployment",
                                              "Bloc",
                                              "Day", 
                                              "Responsetime",
                                              "Device"),
                           Covariate_label = c("Age",
                                              "Woman (vs. man)",
                                              "Lives in capital region",
                                              "Has university degree",
                                              "Annual pers. income more than 400,000 DKK",
                                              "Blue collar worker",
                                              "White collar worker",
                                              "Retired",
                                              "Read news about welfare state at least once a week",
                                              "Fears unemployment",
                                              "Voted for right-wing (vs. left-wing) party",
                                              "Day of interview",
                                              "Survey response time",
                                              "Responded on PC/Mac (vs. Phone/Tablet)"))

Heterogeneity <- Hetereogeneity_results$results %>% 
  left_join(Covariate_lookup, by = c("covar" = "Covariate_name")) %>% 
  mutate(
    Covariate_label = Covariate_label %>%factor() %>% fct_relevel("Age",
                                              "Woman (vs. man)",
                                              "Lives in capital region",
                                              "Has university degree",
                                              "Annual pers. income more than 400,000 DKK",
                                              "Blue collar worker",
                                              "White collar worker",
                                              "Retired",
                                              "Read news about welfare state at least once a week",
                                              "Fears unemployment",
                                              "Voted for right-wing (vs. left-wing) party",
                                              "Day of interview",
                                              "Survey response time",
                                              "Responded on PC/Mac (vs. Phone/Tablet)"),
    Level = Level %>% factor %>% fct_relevel("Non-citizens receive 50% less",
                                             "Non-citizens receive 25% less",
                                              "Lived in DK 10 out of 11 years",
                                              "Lived in DK 9 out of 10 years",
                                              "Lived in DK 3 out of 4 years"))

ggplot(data = Heterogeneity, aes(x = Covariate_label, y = Level, fill = importance)) +
  facet_grid(Attribute ~ ., space = "free", scales = "free", switch = "y") +
  geom_tile() +
  scale_fill_gradient(low="white", high="firebrick2", breaks = c(25,50,75,100)) +
  labs(x = "Respondent Characteristic", y = "Welfare Chauvinistic Policy Element", fill = "Imp.") +
  theme_classic2() +
  theme(text = element_text(size = 14),
        axis.text.x = element_text(angle=45, vjust = 1, hjust = 1))

ggsave("Figure_2.jpg", plot = last_plot(), dpi = 300, width = 10, height = 8)

# Panel B

imce_sched <- IMCE_model$imce

imce_sched$Bloc <- as.factor(imce_sched$Bloc)
attribute <- "Non-citizens receive 25% less"
covar = "Bloc"

plot_data <- imce_sched[!is.na(imce_sched[[covar]]) & !is.na(imce_sched[[attribute]]),]

plot_data <- plot_data[order(plot_data[[attribute]]),]
plot_data$effect_order <- 1:nrow(plot_data)

plot_data$imce <- plot_data[[attribute]]
plot_data$covar <- plot_data[[covar]]

conf_int <- IMCE_model$imce_lower %>% 
  select(Observation, `Non-citizens receive 25% less`) %>% 
  left_join({IMCE_model$imce_upper %>% select(Observation, `Non-citizens receive 25% less`)},
            by = "Observation")

colnames(conf_int) <- c("Observation","lower","upper")

plot_data <- left_join(plot_data, conf_int, by = "Observation")

effect_line <- ggplot(plot_data, aes(x = effect_order, y = imce)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey", alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = mean(plot_data[[attribute]]), 
             color = "dodgerblue",
             size = 0.2) +
  labs(x = "", y = "IMCE") +
  theme_classic2() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

covar_density <- ggplot(plot_data, aes(x = effect_order, fill = covar)) +
  geom_histogram(binwidth = 60,position="stack") +
  labs(x = "", y = "Density", color = str_to_title(covar), fill = "Bloc") +
  scale_fill_manual(values = c("red", "blue")) +
  theme_classic2() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
  theme(legend.position = "bottom") +
  scale_x_continuous(limits = c(0, nrow(plot_data)))

plot_grid(effect_line, covar_density, ncol = 1, align = "v", rel_heights = c(0.6,0.4))
ggsave("Figure_2b.jpg", dpi=400, width = 10, height = 6)
```

# Supplementary Material

\newpage
## Appendix B: Variables
```{r Variables}
Benefits_data %>% 
  select(Age, Gender, Copenhagen, University, Income, Blue_collar, White_collar, Retired, News_welfare, Fear_unemployment, Bloc, Day, Responsetime, Device) %>%
  tbl_summary(
    statistic = list(all_continuous() ~ "{mean} ({sd})",
                     all_categorical() ~ "{p}% {n}")) %>%
  as_kable_extra() %>% 
  save_kable(file = "Table_B1.html")

Benefits_subgroup %>% 
  select(Age, Gender, Copenhagen, University, Income, Blue_collar, White_collar, Retired, News_welfare, Fear_unemployment, Bloc, Day, Responsetime, Device) %>%
  tbl_summary(
    statistic = list(all_continuous() ~ "{mean} ({sd})",
                     all_categorical() ~ "{p}% {n}")) %>%
  as_kable_extra() %>% 
  save_kable(file = "Table_B2.html")
```

\newpage
## Appendix C: Design Diagnostics
### C.1: Balance tables
```{r Balance table, results='asis'}
## Balance test across covariates for attributes of interest
Table_C1a <- Benefits_data %>%
select(Citizenship, Age, Gender, Copenhagen, University, Blue_collar, White_collar, Retired, News_welfare, Fear_unemployment, Day, Responsetime, Device) %>%
tbl_summary(
by = Citizenship,
statistic = list(all_continuous() ~ "{mean} ({sd})",
all_categorical() ~ "{p}% {n}"))
Table_C1a %>%
  as_gt() %>%
gt::gtsave(filename = "Table_C1a.html")


Table_C1b <- Benefits_data %>%
select(Residence, Age, Gender, Copenhagen, University, Blue_collar, White_collar, Retired, News_welfare, Fear_unemployment, Responsetime, Device) %>%
tbl_summary(
by = Residence,
statistic = list(all_continuous() ~ "{mean} ({sd})",
all_categorical() ~ "{p}% {n}")) 
Table_C1b %>%
  as_gt() %>%
gt::gtsave(filename = "Table_C1b.html")
```

\newpage
### C.2: Left vs. right profile choice
```{r left vs. right, results='asis'}
Main_mod_right <- lm_robust(
  data = Benefits_analysis %>% filter(Profile == "Right"), 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod_right)

Main_mod_left <- lm_robust(
  data = Benefits_analysis %>% filter(Profile == "Left"), 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod_left)

htmlreg(
  list(Main_mod_left, Main_mod_right),
  custom.header = list("Main analysis (Left profiles)" = 1, "Main analysis (Right profiles)" =2),
  include.ci = TRUE, include.rmse = FALSE, digits = 3,
  doctype = FALSE, stars = c(0.001, 0.01, 0.05, 0.1), symbol = "+", 
  custom.coef.names = c("(Intercept)", "Must have lived in DK 3 out of 4 years", "Must have lived in DK 9 out of 10 years", "Must have lived in DK 10 out of 11 years", "Non-citizens receive 25% less", "Non-citizens receive 50% less", "Proposed by party from red bloc", "Lower costs than status quo", "Higher costs than status quo", "Disabled receive 25% more", "Disabled receive 50% more", "Parents receive 25% more", "Parents receive 50% more", "No-show at employment agency meeting: 10% sanction for 1 month", "No-show at employment agency meeting: 20% sanction for 1 month", "Must accept job offer that matches skills", "Must attend job training", "Must have worked minimum 225 hours last 12 months","Must have worked minimum 225 hours last 12 months. Rule also applies to spouse", "Must have worked minimum 300 hours last 12 months. Rule also applies to spouse"),
file = "Table_C2.doc", star.symbol = "\\*")

```

\newpage
### C.3: Task-specific AMCEs
```{r Task, results='asis'}
Main_mod_task1 <- lm_robust(
  data = Benefits_analysis %>% filter(Task == "First choice task"), 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod_task1)

Main_mod_task2 <- lm_robust(
  data = Benefits_analysis %>% filter(Task == "Second choice task"), 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod_task2)

Main_mod_task3 <- lm_robust(
  data = Benefits_analysis %>% filter(Task == "Third choice task"), 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod_task3)

Main_mod_task4 <- lm_robust(
  data = Benefits_analysis %>% filter(Task == "Fourth choice task"), 
  formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(Main_mod_task4)

htmlreg(
  list(Main_mod_task1,Main_mod_task2,Main_mod_task3,Main_mod_task4),
  custom.header = list("Main analysis (first task)" = 1, "Main analysis (second task)" =2,"Main analysis (third task)" = 3, "Main analysis (fourth task)" =4),
  include.ci = TRUE, include.rmse = FALSE, digits = 3,
  doctype = FALSE, stars = c(0.001, 0.01, 0.05, 0.1), symbol = "+", 
  custom.coef.names = c("(Intercept)", "Must have lived in DK 3 out of 4 years", "Must have lived in DK 9 out of 10 years", "Must have lived in DK 10 out of 11 years", "Non-citizens receive 25% less", "Non-citizens receive 50% less", "Proposed by party from red bloc", "Lower costs than status quo", "Higher costs than status quo", "Disabled receive 25% more", "Disabled receive 50% more", "Parents receive 25% more", "Parents receive 50% more", "No-show at employment agency meeting: 10% sanction for 1 month", "No-show at employment agency meeting: 20% sanction for 1 month", "Must accept job offer that matches skills", "Must attend job training", "Must have worked minimum 225 hours last 12 months","Must have worked minimum 225 hours last 12 months. Rule also applies to spouse", "Must have worked minimum 300 hours last 12 months. Rule also applies to spouse"),
file = "Table_C3.doc", star.symbol = "\\*")
```

\newpage
### C.4: Row-order effects
```{r Row, results='asis'}
## Testing for potential row order effects
mod_citizenship_row <- lm_robust(
  data = Benefits_analysis, 
  formula = Choice ~Residence + Citizenship*Citizenship_row + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(mod_citizenship_row)

plotdata_C3a <- bind_rows(
  tidy(mod_citizenship_row) %>% filter(
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row2" |   
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row2" |  
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row3" |  
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row3" |  
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row4" | 
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row4" |  
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row5" |  
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row5" |  
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row6" |  
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row6" |  
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row7" |  
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row7" |  
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row8" | 
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row8" | 
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row9" |  
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row9" | 
term == "CitizenshipNon-citizens receive 25% less:Citizenship_row10" |  
term == "CitizenshipNon-citizens receive 50% less:Citizenship_row10"))

plotdata_C3a <- cbind(ID = 1:nrow(plotdata_C3a), plotdata_C3a)


plotdata_C3a <- plotdata_C3a %>%
  mutate(
    min90 = (estimate - qt(0.95, df) * std.error) * 100,
    min95 = (estimate - qt(0.975, df) * std.error) * 100,
    max90 = (estimate + qt(0.95, df) * std.error) * 100,
    max95 = (estimate + qt(0.975, df) * std.error) * 100,
    label_up = round(estimate, digits = 3) * 100,
    estimate = estimate * 100,
    Level = case_when(
    ID == 1 | ID == 3 |ID == 5 |ID == 7 |ID == 9 |ID == 11 |ID == 13 |ID == 15 |ID == 17 ~ "25% cuts for non-citizens",
    ID == 2 | ID == 4 |ID == 6 |ID == 8 |ID == 10 |ID == 12 |ID == 14 |ID == 16 |ID == 18 ~ "50% cuts for non-citizens"),
    Term = case_when(
    ID == 1 | ID == 2 ~ "2 row",
    ID == 3 | ID == 4 ~ "3 row",
    ID == 5 | ID == 6 ~ "4 row",
    ID == 7 | ID == 8 ~ "5 row",
    ID == 9 | ID == 10 ~ "6 row",
    ID == 11 | ID == 12 ~ "7 row",
    ID == 13 | ID == 14 ~ "8 row",
    ID == 15 | ID == 16 ~ "9 row",
    ID == 17 | ID == 18 ~ "10 row") %>% factor() %>% fct_relevel("10 row",
"9 row",
"8 row",
"7 row",
"6 row",
"5 row",
"4 row",
"3 row",
"2 row"))

## Plot   
Figure_C1a <- ggplot(data = plotdata_C3a, aes(y = estimate, x = Term, shape = Level)) + 
  geom_hline(yintercept = 0, color = "#901A1E", lty = "dashed") +
  geom_linerange(aes(ymin = min95, ymax = max95), color = "#808080", position = position_dodge(width = 0.7)) +
  geom_pointrange(aes(ymin = min90, ymax = max90), position = position_dodge(width = 0.7), fill = "white", size = 0.8) +
  geom_text(aes(label = label_up), position = position_dodge(width = 0.7), vjust = -1.2, size = 3) +
  coord_flip() +
  facet_grid(~ Level) +
  labs(x = "", y = "AMCE-difference to row 1") + 
  theme_classic2() +
  theme(legend.position="bottom")

ggsave("Figure_C1a.jpg", plot = Figure_C1a, width = 6, height = 4, dpi = 600)

mod_residence_row <- lm_robust(
  data = Benefits_analysis, 
  formula = Choice ~ Residence*Residence_row + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  clusters = ID,
  weights = Weight)
summary(mod_residence_row)

plotdata_C3b <- bind_rows(
  tidy(mod_residence_row) %>% filter(
term == "ResidenceLived in DK 3 out of 4 years:Residence_row2" |   
term == "ResidenceLived in DK 3 out of 4 years:Residence_row3" |   
term == "ResidenceLived in DK 3 out of 4 years:Residence_row4" |   
term == "ResidenceLived in DK 3 out of 4 years:Residence_row5" |   
term == "ResidenceLived in DK 3 out of 4 years:Residence_row6" |   
term == "ResidenceLived in DK 3 out of 4 years:Residence_row7" |   
term == "ResidenceLived in DK 3 out of 4 years:Residence_row8" |   
term == "ResidenceLived in DK 3 out of 4 years:Residence_row9" |   
term == "ResidenceLived in DK 3 out of 4 years:Residence_row10" | 
term == "ResidenceLived in DK 9 out of 10 years:Residence_row2" |   
term == "ResidenceLived in DK 9 out of 10 years:Residence_row3" |   
term == "ResidenceLived in DK 9 out of 10 years:Residence_row4" |   
term == "ResidenceLived in DK 9 out of 10 years:Residence_row5" |   
term == "ResidenceLived in DK 9 out of 10 years:Residence_row6" |   
term == "ResidenceLived in DK 9 out of 10 years:Residence_row7" |   
term == "ResidenceLived in DK 9 out of 10 years:Residence_row8" |   
term == "ResidenceLived in DK 9 out of 10 years:Residence_row9" |   
term == "ResidenceLived in DK 9 out of 10 years:Residence_row10" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row2" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row3" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row4" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row5" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row6" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row7" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row8" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row9" |   
term == "ResidenceLived in DK 10 out of 11 years:Residence_row10"))
  
plotdata_C3b <- cbind(ID = 1:nrow(plotdata_C3b), plotdata_C3b)


plotdata_C3b <- plotdata_C3b %>%
  mutate(
    min90 = (estimate - qt(0.95, df) * std.error) * 100,
    min95 = (estimate - qt(0.975, df) * std.error) * 100,
    max90 = (estimate + qt(0.95, df) * std.error) * 100,
    max95 = (estimate + qt(0.975, df) * std.error) * 100,
    label_up = round(estimate, digits = 3) * 100,
    estimate = estimate * 100,
    Level = case_when(
    ID == 1 | ID == 4 |ID == 7 |ID == 10 |ID == 13 |ID == 16 |ID == 19 |ID == 22 |ID == 25 ~ "Must have lived 3 out of 4 y in DK",
    ID == 2 | ID == 5 |ID == 8 |ID == 11 |ID == 14 |ID == 17 |ID == 20 |ID == 23 |ID == 26 ~ "Must have lived 9 out of 10 y in DK",
    ID == 3 | ID == 6 |ID == 9 |ID == 12 |ID == 15 |ID == 18 |ID == 21 |ID == 24 |ID == 27 ~ "Must have lived 10 out of 11 y in DK") %>% factor() %>% fct_relevel("Must have lived 3 out of 4 y in DK",
                                                                                                                                                                  "Must have lived 9 out of 10 y in DK",
                                                                                                                                                                  "Must have lived 10 out of 11 y in DK"),
    Term = case_when(
    ID == 1 | ID == 2 | ID == 3 ~ "2 row",
    ID == 4 | ID == 5 | ID == 6 ~ "3 row",
    ID == 7 | ID == 8 | ID == 9 ~ "4 row",
    ID == 10 | ID == 11 | ID == 12 ~ "5 row",
    ID == 13 | ID == 14 | ID == 15 ~ "6 row",
    ID == 16 | ID == 17 | ID == 18 ~ "7 row",
    ID == 19 | ID == 20 | ID == 21 ~ "8 row",
    ID == 22 | ID == 23 | ID == 24 ~ "9 row",
    ID == 25 | ID == 26 | ID == 27 ~ "10 row") %>% factor() %>% fct_relevel("10 row",
"9 row",
"8 row",
"7 row",
"6 row",
"5 row",
"4 row",
"3 row",
"2 row"))

## Plot
Figure_C1b <- ggplot(data = plotdata_C3b, aes(y = estimate, x = Term, shape = Level)) + 
  geom_hline(yintercept = 0, color = "#901A1E", lty = "dashed") +
  geom_linerange(aes(ymin = min95, ymax = max95), color = "#808080", position = position_dodge(width = 0.7)) +
  geom_pointrange(aes(ymin = min90, ymax = max90), position = position_dodge(width = 0.7), fill = "white", size = 0.8) +
  geom_text(aes(label = label_up), position = position_dodge(width = 0.7), vjust = -1.2, size = 3) +
  coord_flip() +
  facet_grid(~ Level) +
  labs(x = "", y = "AMCE-difference to row 1") + 
  theme_classic2() +
  theme(legend.position="bottom")

ggsave("Figure_C1b.jpg", plot = Figure_C1b, width = 8, height = 4, dpi = 600)
```

\newpage
### C.5: AMCEs without speeders and slowpokes
```{r Response time, results='asis'}
# Histogram of questionnaire response time, log scale
Figure_C2 <- ggplot(Benefits_data, aes(x = Responsetime)) +
  geom_histogram(aes(y = stat(count / sum(count))))+
  scale_x_log10() +
  labs(x = "Survey response time (log scale)",
       y = "% responses") +
  theme_classic()

ggsave("Figure_C2.jpg", plot = Figure_C2, width = 6, height = 6, dpi = 400)

# Re-estimate analysis excluding fastest and slowest respondents
Main_mod_responsetime <- lm_robust(
data = Benefits_data %>% filter(Responsetime > 2.5 & Responsetime < 16),
formula = Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
clusters = ID,
weights = Weight)
summary(Main_mod_responsetime)

htmlreg(
list(Main_mod, Main_mod_responsetime),
include.ci = TRUE, include.rmse = FALSE, digits = 3,
doctype = FALSE, stars = c(0.001, 0.01, 0.05, 0.1), symbol = "+", 
custom.header = list("Full sample" = 1, "Without speeders and slowpokes" =2 ),
  custom.coef.names = c("(Intercept)", "Must have lived in DK 3 out of 4 years", "Must have lived in DK 9 out of 10 years", "Must have lived in DK 10 out of 11 years", "Non-citizens receive 25% less", "Non-citizens receive 50% less", "Proposed by party from red bloc", "Lower costs than status quo", "Higher costs than status quo", "Disabled receive 25% more", "Disabled receive 50% more", "Parents receive 25% more", "Parents receive 50% more", "Sanctioning 10% of cash benefits", "Sanctioning 20% of cash benefits", "Need to accept job offer that matches skills", "Need to attend job training", "225 Hours work requirements", "225 hours work requirements for both partners in a union", "300 hours work requirements for both partners in a union"),
caption = "Average marginal component effects: Excluding those who took less than 2.5 and more than 16 minutes",
caption.above = TRUE, scriptsize = TRUE, file = "Table_C4.doc", star.symbol = "\\*")
```

\newpage
### Appendix E: Marginal Means
```{r Marginal means, results='asis'}
Design <- svydesign(
  ids = ~ID,
  weights = ~Weight,
  data = Benefits_analysis
)

Model <- svyglm(
  Choice ~ Residence + Citizenship + Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work, 
  design = Design
)

tidy(Model)

Main_mod_MM <- marginal_means(
  Model, 
  newdata = c(
    "Residence", "Citizenship",
    "Proposer", "Costs", "Handicap", "Parents", "Sanctions", "Joboffer", "Jobtraining", "Work"),
  wts = "cells"
)

## Generate the plotdata
Plotdata_E1 <- tidy(Main_mod_MM) %>% 
  mutate(
    min95 = conf.low * 100,
    max95 = conf.high * 100,
    Attribute = case_when(
      term == "Residence" ~ "Residence",
      term == "Citizenship" ~ "Non-citizens",
      term == "Proposer" ~ "Proposer",
      term == "Costs" ~ "Tax costs",
      term == "Handicap" ~ "Disability",
      term == "Parents" ~ "Parents",
      term == "Sanctions" ~ "Sanctions",
      term == "Joboffer" ~ "Job offer",
      term == "Jobtraining" ~ "Job training",
      term == "Work" ~ "Prior work") %>% factor() %>% fct_relevel("Residence",
                                               "Non-citizens",
                                               "Disability",
                                               "Parents",
                                               "Proposer",
                                               "Tax costs",
                                               "Sanctions",
                                               "Job offer",
                                               "Job training",
                                               "Prior work"),
    label_up = round(estimate, digits = 3) * 100,
    estimate = estimate * 100) %>% drop_na()


Plotdata_E1 <- Plotdata_E1 %>%
  mutate(
    Term = case_when(
      value == "No work requirements" ~ "No work requirements",
      value == "225 hours" ~ "225 h last 12 months",
      value == "225 hours both partners in a union" ~ "225 h last 12 months incl. spouse",
      value == "300 hours both partners in a union" ~ "300 h last 12 months incl. spouse",
      value == "Residency in DK" ~ "Lives in DK",
      value == "Lived in DK 3 out of 4 years" ~ "Lived in DK 3 out of 4 y",
      value == "Lived in DK 9 out of 10 years" ~ "Lived in DK 9 out of 10 y",
      value == "Lived in DK 10 out of 11 years" ~ "Lived in DK 10 out of 11 y",
      value == "Non-citizens receive the same" ~ "Non-citizens receive the same",
      value == "Non-citizens receive 25% less" ~ "Non-citizens receive 25% less",
      value == "Non-citizens receive 50% less" ~ "Non-citizens receive 50% less",
      value == "Party from blue bloc" ~ "Right-wing party",
      value == "Party from red bloc" ~ "Left-wing party",
      value == "Lower personal costs than previous reform" ~ "Pay 0.5% less in annual taxes",
      value == "Same personal costs as previous reform" ~ "No change in annual taxes",
      value == "Higher personal costs than previous reform" ~ "Pay 0.5% more in annual taxes",
      value == "Disabled receive the same" ~ "Disabled receive the same",
      value == "Disabled receive 25% more" ~ "Disabled receive 25% more",
      value == "Disabled receive 50% more" ~ "Disabled receive 50% more",
      value == "Parents receive the same" ~ "Parents receive the same",
      value == "Parents receive 25% more" ~ "Parents receive 25% more",
      value == "Parents receive 50% more" ~ "Parents receive 50% more",
      value == "No sanctions" ~ "No sanctions",
      value == "Receive 10% less the following month" ~ "Receive 10% less for 1 month",
      value == "Receive 20% less the following month" ~ "Receive 20% less for 1 month",
      value == "Not required to accept joboffer that matches skills" ~ "No job offer requirements",
      value == "Required to accept joboffer that matches skills" ~ "Must accept job offer",
      value ==  "Not required to participate in job training" ~  "No job training requirements",
      value == "Required to participate in job training" ~ "Must attend job training") %>% factor() %>% fct_relevel("300 h last 12 months incl. spouse",
                                                                                              "225 h last 12 months incl. spouse",
                                                                                              "225 h last 12 months",
                                                                                              "No work requirements",
                                                                                              "Must attend job training",
                                                                                              "No job training requirements",
                                                                                              "Must accept job offer",
                                                                                              "No job offer requirements",
                                                                                              "Receive 20% less for 1 month",
                                                                                              "Receive 10% less for 1 month",
                                                                                              "No sanctions",
                                                                                              "Left-wing party",
                                                                                              "Right-wing party",
                                                                                              "Pay 0.5% more in annual taxes",
                                                                                              "No change in annual taxes",
                                                                                              "Pay 0.5% less in annual taxes",
                                                                                              "Parents receive 50% more",
                                                                                              "Parents receive 25% more",
                                                                                              "Parents receive the same",                                                                                              
                                                                                              "Disabled receive 50% more",
                                                                                              "Disabled receive 25% more",
                                                                                              "Disabled receive the same",                                                                                               
                                                                                              "Non-citizens receive 50% less",
                                                                                              "Non-citizens receive 25% less",
                                                                                              "Non-citizens receive the same",
                                                                                              "Lived in DK 10 out of 11 y",
                                                                                              "Lived in DK 9 out of 10 y",
                                                                                              "Lived in DK 3 out of 4 y",
                                                                                              "Lives in DK")) %>% select(Attribute, Term, estimate, label_up, min95, max95)


## Plot   
Figure_E1 <- ggplot(data = Plotdata_E1, 
       aes(y = estimate, x = Term)) +
  geom_hline(yintercept = 50, color = "#901A1E", lty = "dashed") +
  geom_errorbar(aes(ymin = min95, ymax = max95), width = 0.2, color = "#000000") +  # Use geom_errorbar
  geom_point(size = 0.5) +
  geom_text(aes(label = label_up), vjust = -1.2, hjust = -0.5, size = 2.5) +
  scale_y_continuous(breaks = c(40,45,50,55,60),
                      labels = c("40", "45", "50", "55", "60")) +
  coord_flip() +
  labs(x = "", 
       y = "Marginal Means") + 
  facet_grid(Attribute~., scales = "free_y", switch = "y") +
  theme_classic2()


## Save Figure_E1
ggsave("Figure_E1.jpg", plot = last_plot(), width = 8, height = 10, dpi = 600)
```

\newpage
## Appendix F: Sub-group effects
### F.1 Treatment-by-covariate interactions 
```{r Sub-group effects - LOR x Gender, results='asis'}
imce_sched$Gender <- as.factor(imce_sched$Gender)
attribute <- "Lived in DK 3 out of 4 years"
covar = "Gender"

plot_data <- imce_sched[!is.na(imce_sched[[covar]]) & !is.na(imce_sched[[attribute]]),]

plot_data <- plot_data[order(plot_data[[attribute]]),]
plot_data$effect_order <- 1:nrow(plot_data)

plot_data$imce <- plot_data[[attribute]]
plot_data$covar <- plot_data[[covar]]

conf_int <- IMCE_model$imce_lower %>% 
  select(Observation, `Lived in DK 3 out of 4 years`) %>% 
  left_join({IMCE_model$imce_upper %>% select(Observation, `Lived in DK 3 out of 4 years`)},
            by = "Observation")

colnames(conf_int) <- c("Observation","lower","upper")

plot_data <- left_join(plot_data, conf_int, by = "Observation")

effect_line <- ggplot(plot_data, aes(x = effect_order, y = imce)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey", alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = mean(plot_data[[attribute]]), 
             color = "dodgerblue",
             size = 0.2) +
  labs(x = "", y = "IMCE") +
  theme_classic2() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

covar_density <- ggplot(plot_data, aes(x = effect_order, fill = covar)) +
  geom_histogram(binwidth = 60,position="stack") +
  labs(x = "", y = "Density", color = str_to_title(covar), fill = "Gender") +
  scale_fill_manual(values = c("blue", "pink")) +
  theme_classic2() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
  theme(legend.position = "bottom") +
  scale_x_continuous(limits = c(0, nrow(plot_data)))

plot_grid(effect_line, covar_density, ncol = 1, align = "v", rel_heights = c(0.6,0.4))
ggsave("Figure_F1.jpg", dpi=400, width = 10, height = 6)

attribute <- "Lived in DK 9 out of 10 years"
covar = "Gender"

plot_data <- imce_sched[!is.na(imce_sched[[covar]]) & !is.na(imce_sched[[attribute]]),]

plot_data <- plot_data[order(plot_data[[attribute]]),]
plot_data$effect_order <- 1:nrow(plot_data)

plot_data$imce <- plot_data[[attribute]]
plot_data$covar <- plot_data[[covar]]

conf_int <- IMCE_model$imce_lower %>% 
  select(Observation, `Lived in DK 9 out of 10 years`) %>% 
  left_join({IMCE_model$imce_upper %>% select(Observation, `Lived in DK 9 out of 10 years`)},
            by = "Observation")

colnames(conf_int) <- c("Observation","lower","upper")

plot_data <- left_join(plot_data, conf_int, by = "Observation")

effect_line <- ggplot(plot_data, aes(x = effect_order, y = imce)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey", alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = mean(plot_data[[attribute]]), 
             color = "dodgerblue",
             size = 0.2) +
  labs(x = "", y = "IMCE") +
  theme_classic2() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

covar_density <- ggplot(plot_data, aes(x = effect_order, fill = covar)) +
  geom_histogram(binwidth = 60,position="stack") +
  labs(x = "", y = "Density", color = str_to_title(covar), fill = "Gender") +
  scale_fill_manual(values = c("blue", "pink")) +
  theme_classic2() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
  theme(legend.position = "bottom") +
  scale_x_continuous(limits = c(0, nrow(plot_data)))

plot_grid(effect_line, covar_density, ncol = 1, align = "v", rel_heights = c(0.6,0.4))
ggsave("Figure_F2.jpg", dpi=400, width = 10, height = 6)

attribute <- "Lived in DK 10 out of 11 years"
covar = "Gender"

plot_data <- imce_sched[!is.na(imce_sched[[covar]]) & !is.na(imce_sched[[attribute]]),]

plot_data <- plot_data[order(plot_data[[attribute]]),]
plot_data$effect_order <- 1:nrow(plot_data)

plot_data$imce <- plot_data[[attribute]]
plot_data$covar <- plot_data[[covar]]

conf_int <- IMCE_model$imce_lower %>% 
  select(Observation, `Lived in DK 10 out of 11 years`) %>% 
  left_join({IMCE_model$imce_upper %>% select(Observation, `Lived in DK 10 out of 11 years`)},
            by = "Observation")

colnames(conf_int) <- c("Observation","lower","upper")

plot_data <- left_join(plot_data, conf_int, by = "Observation")

effect_line <- ggplot(plot_data, aes(x = effect_order, y = imce)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey", alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_hline(yintercept = mean(plot_data[[attribute]]), 
             color = "dodgerblue",
             size = 0.2) +
  labs(x = "", y = "IMCE") +
  theme_classic2() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

covar_density <- ggplot(plot_data, aes(x = effect_order, fill = covar)) +
  geom_histogram(binwidth = 60,position="stack") +
  labs(x = "", y = "Density", color = str_to_title(covar), fill = "Gender") +
  scale_fill_manual(values = c("blue", "pink")) +
  theme_classic2() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  guides(fill=guide_legend(nrow=1,byrow=TRUE)) +
  theme(legend.position = "bottom") +
  scale_x_continuous(limits = c(0, nrow(plot_data)))

plot_grid(effect_line, covar_density, ncol = 1, align = "v", rel_heights = c(0.6,0.4))
ggsave("Figure_F3.jpg", dpi=400, width = 10, height = 6)
```

\newpage
### F.2 Treatment-by-treatment interactions
```{r TxT interactions}
# Define the interaction terms
interaction_terms <- c("Proposer", "Costs", "Handicap", "Parents", "Sanctions", "Joboffer", "Jobtraining", "Work")

# Function to run lm_robust models for a given interaction variable
run_models <- function(interaction_var, primary_var) {
  formula_str <- paste("Choice ~", primary_var, "+", interaction_var, "*", primary_var, "+ Proposer + Costs + Handicap + Parents + Sanctions + Joboffer + Jobtraining + Work")
  model <- lm_robust(as.formula(formula_str), data = Benefits_analysis, clusters = ID, weights = Weight)
  
  # Store results in a tidy format
  tidy(model) %>% 
    mutate(interaction = interaction_var, primary = primary_var)
}

# Run models for both "Citizenship" and "Residence" interactions
results_citizenship <- map_dfr(interaction_terms, run_models, primary_var = "Citizenship")
results_residence <- map_dfr(interaction_terms, run_models, primary_var = "Residence")

# Combine results into a single dataframe and filter only interaction terms
Interaction_results <- bind_rows(results_citizenship, results_residence) %>% 
  filter(str_detect(term, ":"))

Interaction_results <- Interaction_results %>%
  mutate(
    Estimate = estimate*100,
    Min95 = conf.low*100,
    Max95 = conf.high*100,
    Term = case_when(
      grepl("ProposerParty from red bloc", term) ~ "Proposer is left-wing party",
      grepl("CostsLower personal costs than previous reform", term) ~ "Pay 0.5% less in annual taxes",
      grepl("CostsHigher personal costs than previous reform", term) ~ "Pay 0.5% more in annual taxes",
      grepl("HandicapDisabled receive 25% more", term) ~ "Disabled receive 25% more",
      grepl("HandicapDisabled receive 50% more", term) ~ "Disabled receive 50% more",
      grepl("ParentsParents receive 25% more", term) ~ "Claimants with at-home kids receive 25% more",
      grepl("ParentsParents receive 50% more", term) ~ "Claimants with at-home kids receive 50% more",
      grepl("SanctionsReceive 10% less the following month", term) ~ "Receive 10% less for 1 month",
      grepl("SanctionsReceive 20% less the following month", term) ~ "Receive 20% less for 1 month",
      grepl("JobofferRequired to accept joboffer that matches skills", term) ~ "Must accept job offer",
      grepl("JobtrainingRequired to participate in job training", term) ~ "Must attend job training program",
      grepl("Work225 hours$", term) ~ "225 h last 12 months",
      grepl("Work225 hours both partners in a union", term) ~ "225 h last 12 months incl. spouse",
      grepl("Work300 hours both partners in a union", term) ~ "300 h last 12 months incl. spouse"),
    Wc = case_when(
      grepl("CitizenshipNon-citizens receive 25% less", term) ~ "Non-citizens receive 25% less",
      grepl("CitizenshipNon-citizens receive 50% less", term) ~ "Non-citizens receive 50% less",
      grepl("ResidenceLived in DK 3 out of 4 years", term) ~ "Lived in DK 3 of 4 years",
      grepl("ResidenceLived in DK 9 out of 10 years", term) ~ "Lived in DK 9 of 10 years",
      grepl("ResidenceLived in DK 10 out of 11 years", term) ~ "Lived in DK 10 of 11 years") %>% factor() %>%fct_relevel("Non-citizens receive 25% less",
                                                                                                                          "Non-citizens receive 50% less",
                                                                                                                          "Lived in DK 3 of 4 years",
                                                                                                                          "Lived in DK 9 of 10 years",
                                                                                                                          "Lived in DK 10 of 11 years")) %>% select(Estimate, Min95, Max95, Term, Wc, interaction, primary)

Figure_interaction <- ggplot(Interaction_results, aes(x = Estimate, y = Term)) +
  geom_point(size = 3) +  # Effect size points
  geom_errorbarh(aes(xmin = Min95, xmax = Max95), height = 0.3) +  # Confidence intervals
  geom_vline(xintercept = 0, color = "#901A1E", lty = "dashed") +
  theme_classic2() +
  coord_flip() +
  labs(x = "Extra Effect of Other Attributes",
       y = "") +
  facet_wrap(~Wc, scales = "free_y") +
  theme(axis.text.y = element_text(size = 8),
        axis.text.x = element_text(size = 8, angle=70, vjust = 1, hjust = 1))


ggsave("Figure_F4.jpg", plot = last_plot(), width = 14, height = 12, dpi = 600)
```
